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Abstract 

In the analysis of composite materials with heterogeneous microstructures, full res- 
olution of the heterogeneities using classical numerical approaches can be computa- 
tionally prohibitive. This paper presents a micromechanics-enhanced finite element 
formulation that accurately captures the mechanical behaviour of heterogeneous ma- 
terials in a computationally efficient manner. The strategy exploits analytical solu- 
tions derived by Eshelby for ellipsoidal inclusions in order to determine the mechan- 
ical perturbation fields as a result of the underlying heterogeneities. Approximation 
functions for these perturbation fields are then incorporated into a finite element 
formulation to augment those of the macroscopic fields. A significant feature of 
this approach is that the finite element mesh does not explicitly resolve the hetero- 
geneities and that no additional degrees of freedom are introduced. In this paper, 
hybrid- Trefftz stress finite elements are utilised and performance of the proposed for- 
mulation is demonstrated with numerical examples. The method is restricted here 
to elastic particulate composites with ellipsoidal inclusions but it has been designed 
to be extensible to a wider class of materials comprising arbitrary shaped inclusions. 
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1. Introduction 

In the analysis of materials with complex microstructures, full resolution of 
the heterogeneities using classical numerical approaches such as the Finite Element 
method can be computationally prohibitive. To overcome this, one option is to model 
the macroscale problem using equivalent properties; however, this can lead to a criti- 
cal loss of information about the finer scale behaviour and poor understanding of the 
heterogeneities' infiuence on the macroscale response. Numerical approaches such 
as computational homogenization (often called FE^) provide an alternative strat- 
egy [HE]. These techniques comprise nested Finite Element analyses, where each 
macroscopic material point response is determined via the numerical solution of an 
RVE subject to the macroscopic strains. Although such approaches have significant 
potential for certain classes of problems, they are still computationally demanding 
and are restricted to situations involving clear separation of scales. 

The objective of this work is to develop a Finite Element formulation for modelling 
the macroscopic mechanical problem that is enhanced to capture the infiuence of the 
underlying heterogeneities. In our approach, the Finite Element mesh is not required 
to explicitly resolve the heterogeneities. Closed-form expressions for the perturbation 
of the mechanical fields due to the presence of the heterogeneities are determined 
and these are then utilised to enhance the Finite Element formulation. 

The ability to capture the eflFect of microstructural features independently of the 
underlying finite element mesh has been an ongoing challenge in computational me- 
chanics research. Partition of Unity methods jSj HI EJ [Gj provide a potential solution 
to this problem, without mesh refinement, by extending a given solution space with 
additional functions and has been successfully applied to problems such as cracks 
and material interfaces. The application of this approach in the context of the cur- 
rent work will be briefiy discussed in this paper, whereby the closed-form solutions 
derived for the mechanical perturbation fields are used to extend the classical finite 
element method. However, it will be shown that there are some disadvantages to this 
approach for the particular problem at hand and an alternative approach, centred 
on the Hybrid- Trefftz stress element formulation [7j, represents the main focus of 
this paper. This method does not result in additional degrees of freedom, although 
it does involve an additional, albeit relatively minor, computational overhead. 

The heterogeneities, although currently restricted to simple shapes (ellipsoids), 
can be randomly sized and randomly distributed without reference to the finite el- 
ement mesh. Therefore, the proposed approach has the potential to be applica- 
ble to a wide range of composite materials, such as fibre reinforced composites [8j, 
porous media [9l [10], functionally graded materials [llj, etc. Moreover, it can be 
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extended to general inclusion shapes by evaluating the perturbation functions nu- 
merically p[2l[T3]. 

The paper is structured as follows. The methodology of the proposed strategy is 
described in Section [2} Construction of the perturbation approximation functions 
for Finite Element Analysis is derived in Section [3} The implementation into the 
Hybrid-Trefftz stress element formulation containing an arbitrary number of inclu- 
sions is presented in Section [4} Section [5] comprises examples demonstrating the 
model's performance. Finally we present the conclusions as well as a discussion on 
future research directions. An appendix is included that highlights some important, 
but rather technical, aspects of the proposed technique in order to keep the paper 
self-contained. 

2. Micromechanics approach 

This section outlines the strategy to calculate the perturbation of mechanical 
fields due to a heterogeneous microstructure, exploiting the Equivalent Inclusion 
Method [11] in conjunction with analytical micromechanics. Our objective is to 
convert the heterogeneous problem into an equivalent homogeneous problem and to 
derive analytical expressions for the perturbations of the stress, strain and displace- 
ment fields that we can then utilise within a finite element formulation. 

Consider a body consisting of clearly distinguishable heterogeneities in a matrix 
(Fig. nk) subjected to a displacement u and traction t field. The stiffness of such a 




(a) (b) 

Figure 1: Principle of Equivalent Inclusion Method : a) composite body with inclusions, b) homo- 
geneous reference body with additional equivalent eigenstrains 

material is decomposed as follows [14[ [T3] 

c = c° + yc* (1) 
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where C° is the stiffness tensor of a homogeneous matrix and C* = X]f^[C^ — C°] 
is due to the presence of inclusions. C* is nonzero only within the domain Q = 
n^U"'Un^,so that 



As a result of the heterogeneities, the mechanical fields (displacement, strains, stresses) 
experience a perturbation for which we will derive closed-form expressions based on 
analytical micromechanics. Symbolically, we can express the decomposition of the 
mechanical fields as follows: 



where, the superscript indicates the macroscopic component of the fields in the ab- 
sence of heterogeneities and superscript •* indicates the perturbation (or microscopic ) 
component due to the presence of the heterogeneities. It is worth noting that, tra- 
ditionally, in analytical micromechanics, the macroscopic fields are assumed to be 
uniform across the domain, e.g. jl5l[l6]. Here it is assumed that they can be position 
dependent functions of the Neumann and Dirichlet boundary conditions. 

The perturbation fields are determined by employing the equivalent inclusion 
method for a single heterogeneity embedded in a matrix and then extended here 
for multiple heterogeneities. In the equivalent inclusion method, the heterogeneous 
solid is replaced by an equivalent homogeneous solid with uniform material stiffness 

everywhere (Fig. [l^, b) and suitable stress-free eigenstrains ej applied in the 
inclusions so that the homogeneous equivalent solid has the same mechanical 
fields as the original heterogeneous solid. 

2.1. Equivalent inclusion method for single heterogeneity problem 

Consider first a single heterogeneity embedded in a matrix. Following Eshelby's 
fundamental work , this problem can be decomposed into two problems of known 
solution and then assembled back via superposition [HI [17], see Fig. [2| In brief, 
the solution of the inhomogeneity problem requires the determination of the trans- 
formation eigenstrain that induces the identical local mechanical response as the 
original heterogeneous body. Due to the absence of other inclusions, no strain or 
stress concentrations are induced and thus remains constant [17j. 

In the original heterogeneous body, the stress state can be expressed as 




(2) 



(3) 



<T = (T° + <T* = C : [£° + £*] 



(4) 
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Figure 2: Equivalent Inclusion Method: a) inhomogeneity problem, b) infinite homogeneous body, 
c) homogeneous inclusion problem 



In the homogeneous solid, we add a stress- free eigenstrain inside the domain of 
the inclusion, which has the same material stiffness C° as the matrix such that 

o- = C° : + s* - s"] (5) 

It should be noted that s'^ = in the matrix. 

Given that the macroscopic stress is = C° : s^, it can be see from equat- 
ing Eqs. Q and ([5]), that the stress perturbation in the homogeneous solid can be 
expressed as 

0-* = C° : [s* - s"] (6) 
Furthermore, equating Eqs. Q and ^ also results in the following expression: 

C : [s V s*] = C° : + s* - e^] (7) 

where the transformation eigenstrain is as yet unknown. Eshelby's solution of the 
homogeneous inclusion problem [17j, relates the eigenstrain to the perturbation strain 
as follows 

s* = S : (8) 

where S denotes the Eshelby tensor and is a function of the heterogeneity's geometry 
and the material stiffness of both the matrix and heterogeneity. Substituting this 
expression into Eq. ([7| and rearranging yields 

[C - C^] : s° = [C^ : S - C : S - C°] : (9) 

This can be recast in a compact form to give an expression for the eigenstrain 
which depends on the homogeneous strain s^, the material stiffness of both the matrix 
and heterogeneity and the Eshelby tensor as 

= B : s° (10) 
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where the tensor B is provided by: 



B = - [C* : S + C°] ^ : C* (11) 

Once the transformation eigenstrain has been determined, the stress perturbation 
can be computed from Eq. ([g]) in the form 

0-* = C° : [S - I] : B : s° (12) 

It can be seen that the stress perturbation depends on stiffness of the different 
material phases, the macroscopic strain field and the geometry of the heterogeneity. 
This closed-form expression for the stress perturbation is at the heart of the proposed 
finite element enrichment to be discussed later in this paper. It is also useful to derive 
an expression for the displacement perturbation field as follows 

u* = £ : = £ : B : s° (13) 

where the operator C is Sl third order tensor, mapping s'^ ^ u*. For the sake 
of conciseness and to keep the paper self-contained, the detailed derivation of this 



operator can be found in Appendix Appendix A 



2.2. Self- compatibility algorithm for multiple inclusions 

In the case of multiple inclusions, the mechanical perturbation fields within in- 
dividual inclusion domains are no longer uniformly distributed as a result of their 
mutual interaction. Here, we account for this approximately by assuming that the 
eigenfields to be piecewice constant within each inclusion. Thus, the perturbations 
are determined from the Eshelby solution for each individual inclusion, as described 
above, plus an iterative self-compatibility procedure (Tab. [T]) to ensure that the 
solution correctly refiects the infiuence of multiple heterogeneities. 

This procedure iteratively enforces compatibility (see eg. ^18j for further refer- 
ence) between the imposed macroscopic strain and the average eigenstrain inside 
any given inclusion so as to account for the infiuence of the remaining inclusions 
N\i. An iterative algorithm has been chosen because a closed form solution for 
multiple inclusions does not exist and a numerical solution would be prohibitively 
expensive [18j. 



First, the eigenstrain ej for each inclusion i is calculated (Eq. 10) without refer- 
ence to the other inclusions (Line 2). Next, the associated perturbation strain s* 
for each inclusion i is evaluated (Eq. [s]) at the centre of all other inclusions (Line 
3). The mutual interaction of inclusions is then taken into account via a correction 
of the eigenstrain (As[). For each inclusion i, this correction is calculated (Line 8) 
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from the inverse of the inclusion's Eshelby tensor S^~^ and the perturbation strains 
of aU other inclusions, evaluated at the centroid of inclusion i. The perturbation 
strain resulting from inclusion j at the centre of inclusion i is denoted as s*^-. This is 
demonstrated in Fig. |3] for a two inclusion problem in ID. The eigenstrain correction 
As[ is then used to calculate the correction to the perturbation strains (Line 10). 
The algorithm continues until a small Euclidean norm between the last two itera- 
tions of the total eigenstrains is achieved. At convergence, the corresponding stress 
and displacement perturbations are recalculated from the corrected transformation 
eigenstrains. 

It is worthwhile noting that the algorithm does not depend on a particular se- 
quence of inclusions, as follows from the elastic reciprocity theorem fT8| and refer- 
ences therein]. Moreover, since the perturbation fields are calculated for the entire 
macroscopic domain (no RVE is considered), stress admissibility and strain compat- 
ibility, in the sense of macro-micro field relations, are fulfilled a priori. 

The computational complexity of the Self-compatibility algorithm is 0{N'^). How- 
ever, this can be improved by taking into account only those inclusions which have 
a non-negligible infiuence to the inclusion of interest i. Preliminary studies have 
shown this to give a significant computational speed-up and will be reported in a 
future paper. 
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Self Compatibility Algorithm (s^, B^, S^, S^" 

For {i<N) 

< = B. : 

< = : ej 
Set As* = s* 

EndFor 
Do 

For {i < N) 



N) 



N 



As 



s- 

ej = ej + AeJ 
As* = S, : As[ 
EndFor 

While (Ef l|A<ll >^) 



(Eq. 



8) 



Table 1: Self-compatibility algorithm. Note, that r] stands for an acceptable tolerance. 
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Figure 3: Principle of self-compatibility algorithm for double inclusion problem in ID. 

3. Construction of perturbation approximation functions for FEA 

The above methodology can be utilised to formulate an enhanced Finite Element 
formulation. The primary task is to determine appropriate approximation functions 
for the mechanical perturbation fields u*, s* and cr* based on the analytical microme- 
chanics developed above and which can then augment the standard macroscopic field 
approximations. It should be noted that the Voigt-Mandel notation is exclusively 
used in the forthcoming text. 

The perturbation field approximation functions are determined a priori as a linear 
combination of the perturbation fields evaluated analytically for six load cases, with 
self-equilibrium enforced by means of the self-compatibility algorithm outlined above 
(Tab.[l|). Each load case corresponds to a unit component of the macroscopic strain 
vector 6^, i = 1, ... ,6 and the resulting analytically determined stress, strain and 
displacement perturbation fields are arranged, column- by-column into ^g^e? ^6x6 
Ugxg matrices, respectively: 

5* = [ V* ... V* ] , e* = [ ^e* ... ^e* ] 
u* = [ ^u* ... ^u* ] 
where the left superscript refers to a specific load case, 1 to 6. 

3.1. Partition of unity method 

Partition of Unity (PU) Methods (for example the eXtended Finite Element 
Method) extend the underlying basis functions used for interpolating the displace- 
ment field by adding an appropriate set of additional functions. Following |19[ [20] it 
has been shown that the displacement field u(x) within an element can be interpo- 
lated by 

n 

u(x) = 5]](N^(x)a' + N'N^(x)b') (16) 

i=l 
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where n is the number of nodes per element, N^(x) = A^^(x)l is the standard matrix 
of element shape functions for node I is the identity matrix and the standard 
displacement degrees of freedom at node i. is a matrix containing the additional 
basis terms and are the associated additional degrees of freedom at node i. It is 
important to recognise that the element shape functions form a partition of unity, 
i.e. 

EiV^(x) = l (17) 

i=l 

The six analytically derived displacement perturbation functions contained in u* can 
be used as the additional functions to augment the standard basis functions. Thus 



N,(x) 



^Mt(x) ... V(x) ... ... 

... i«^(x) ... S^(x) ... 
... ... ^ul{x) ... S^(x) 



(18) 



With this at hand, the PU-based finite element formulation can be derived, see for 
example pTj. PU methods are particularly useful in problems where the extension of 
the basis functions is introduced on a node by node basis, so that additional degrees 
of freedom are only introduced at nodes where the basis is extended. One obvious 
example of such a local feature that can be modelled in this way is discrete cracks [21 J . 
However, this favourable property is not exploited here because we wish to model a 
large number of heterogeneities throughout the domain. In 3D problems, there are 
3 standard displacement degrees of freedom per node; this would be extended by an 
additional 18 degrees of freedom per node and per heterogeneity with the proposed 
approach. 

It is also worth noting that for standard finite elements, the volume integration 
of the discrete system of equations is relatively straightforward. However, extension 



of the basis functions to include the perturbation functions in Eq. (18) makes this 
process significantly more arduous. For these reasons, an alternative Finite Element 
approach using Hybrid Trefftz Stress elements |22l[7j is considered where the standard 
basis function is not extended, as with PU methods, but enhanced such that no 
additional degrees of freedom are introduced. This is described in the next section. 

4. Hybrid- Trefftz stress element formulation 

In this section a finite element formulation based on an enhancement of a hybrid- 
TreflFtz stress (HTS) element formulation \J] is presented. 
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Figure 4: Elastic body representing HTS element 



The problem requires a solution to the displacement u and stress a fields as a 
result of given boundary displacements u and tractions q on and F^, respectively. 
The displacement and stress fields must fulfil the following governing equations: 



(19) 



where a and e are the column matrix representation of the second order stress and 
strain tensor, respectively, u represents the displacement vector, C is the matrix 
representation of fourth order stiffness tensor and finally U and t represent the applied 
displacements and tractions, respectively. The gradient operator L and the matrix 
of directional cosines N of the outward normal to element boundary F^ have the 
following forms [23] 
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4.1. Stress, strain and displacement approximations 

The macroscopic stress field within the HTS element is approximated as 

(j° = S°v in 



(21) 
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where v is the vector of generahsed stress degrees of freedom, S° denotes the matrix of 
stress approximation functions chosen so as to automaticaUy satisfy the equihbrium 



conditions Eq. (19)^'*^. Thus 



and 



NS?.v on 



(22) 



(23) 



where t represents the traction vector induced by the macroscopic stress approxima- 
tion field. 

The macroscopic strain and displacement fields are expressed analogously to 



Eq. (21) as 



e'' = Ey and u° = U°v (24) 
where and are directly associated with the stress approximation by means of 



the compatibility equation (|19)^ and constitutive equation (19)^ as 

C°LU° 



S?, = C°E° 



(25) 



Since the stress approximation functions are typically polynomial functions, the 
integration of to get is relatively straightforward. 

Rather than extend the solution space to capture the infiuence of the hetero- 



geneities, as was briefiy described in Section |3.1[ here we enhance the macroscopic 
approximations to include the infiuence of the heterogeneities, thereby not increas- 
ing the number of unknowns. The total stress (macroscopic plus perturbation) field 
within the HTS element is approximated, following Eq. ([3]), as 



a = a 



a* = (SS + S:) V in 



(26) 



where S* is the perturbation counterpart to S°. S* can be constructed from Eq. (14): 



a 



(27) 



where 5* is the set of analytically defined stress perturbations for six load cases. 



each one representing a unit component of the macroscopic strain Eq. (14). For 



the purposes of constructing S*, we approximate the macroscopic strain field as 
constant within each finite element and computed as the volume average of the 



actual macroscopic strain field. From Eq. (24), 



^0,ave 



1 



1^^ 



1 



(28) 
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Substituting e^'^^^ for into Eq. (27) leads to 

a* =5*Ef^v (29) 

Thus, the matrix of stress perturbation approximation functions is: 

S:=5*Er (30) 

Analogously, the approximation of total strain and displacement fields within the 
element domain is given by 

e=(E° + E:)v and u = (U° + U;) v (31) 

where the perturbation approximation matrices U* and E* are , as with their macro- 
scopic counterparts, directly associated with the stress approximation as 

S: = C°E; = C°LU: (32) 

It is worthwhile noting that the stress perturbation fields and the corresponding 
traction perturbation fields, approximated as a* = S*v and t* = NS*v, remain in 
self-equilibrium. 

4-2. Static boundary conditions 



Contrary to general condition in Eq. ( 19) , the equilibrium on the element traction 



boundary is imposed only in the weighted residual sense as: 

/ \Nj [N + a*) - t] dP = (33) 
along with Wi representing the matrix of weighting functions. Replacing the total 



stress field in Eq. (33) by its approximation from Eq. (26), the traction boundary 
condition becomes 

/ W[N (S^ + S;) vdP = / W^tdP (34) 

4-3. Kinematic boundary conditions 

Compatibility inside the element domain is also enforced in a weighted residual 
sense, such that: 

/ Wj (e° + e* - Lu) dQ' = (35) 
12 



Next, utilising integration by parts and applying Green's theorem to WjLu, Eq. (35) 
results in 



- [ (NW2)^udP= / (NW2)^udP (36) 
With the current formulation, it is not necessarily possible to find a solution to both 



Eqs. (34) & (36) that satisfies both traction and kinematic boundary conditions act- 



ing on the element boundary. As a consequence, Eq. (36) is relaxed by introducing an 



additional and independent approximation of displacements on the element traction 
boundary: 

ur = Urq in (37) 

Here, q stands for the set of displacement unknowns and Ur is the matrix of boundary 
displacement approximation functions. Such a formulation of the stress element leads 
to a hybrid approach O [22] . 



Given the above consideration, introducing Eq. (37) into Eq. (36) results in 



Wj (E° + E:)vdr2^ 



(L"^W2)"^u dfl" 



I (NW2)^urdP= I (NW2)"^udP 

Jvt Jvf, 



(38) 



4^4' Weighting functions 

In order to achieve an energy-consistent formulation, it is required that all weighted 
terms within the integrals defined above have the dimension of work. The weighting 



functions then directly follow from the integrands in Eq. (33) and Eq. (35) represent- 



ing the increment of internal work of strains within and external work of tractions 
on r^, respectively. These functions thus admit the following forms: 



Wi = Ur, W2 = s: 



(39) 



First, introducing Eq. (39)^ into Eq. (38) and taking into account condition ( [22 
yields 



/ (S°)^(E° + E:)vdn^- / (NS°)^Urqdr= / (NS°)^ndP 



(40) 
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Second, introducing Eq. (39)^ into the traction boundary condition (34) yields 



UfN(SS + S:)vdP = / UftdP 



(41) 



Combining Eqs. (40) & (41) results in a coupled system of linear equations that can 
be expressed in compact form as 

F -A^ 1 r V 1 _ r Pu 1 , . 

-(A + A*) q -p. ^^'^ 



where the submatrices on the left-hand side foUow from Eqs (40, 41) and are given 
by the following integrals 



F 
A 



= / (S°)^(E° + E:) dQ<== / N(S°)^(U° + U:) dP (43) 
= [ UrNS°dP and A* = / UrNS*dP (44) 

and for the terms on the right-hand side it holds 

p, = / (NSS)^ndP, and p,= / UjtdP (45) 

Note that Eq. (43) illustrates that the F matrix can be evaluated via a boundary 
rather than volume integral. Thus all terms in Eq. (42) can be evaluated using 
boundary integrals only. 

The size of the system of equations to be solved simultaneously can be reduced 
via static condensation, representing a significant reduction in computational eflFort. 
First, from the first equation in Eq. (42), the generalised stress degrees of freedom v 
are expressed in terms of the displacement degrees of freedom q as 

v = F-HPu + ATq) (46) 

This is then substituted into the second equation of Eq. (42) to yield: 

[(A + A*)F-1aT] q = - (A + A*)F-Vu (47) 

This sparse system of equations is then solved for the displacement degrees of freedom 
q. Subsequently, the stress degrees of freedom v can then be calculated on an element- 
by-element basis. 

Our implementation of these HTS elements for composite materials (C-HTS el- 
ements) utilises displacement degrees of freedom that are associated with element 
faces rather than vertices. This has the advantage that the bandwidth of the stiffness 
matrix is minimised, as is interprocessor communication. 
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5. Numerical Examples 

The performance of the key components of the proposed strategy (micromechan- 
ical solution, self-compatibihty algorithm, finite element analysis convergence, etc.), 
in terms of efficiency and accuracy, have been explored through two numerical ex- 
amples. 

5.1. Three ellipsoidal inclusions in matrix 

This example comprises three ellipsoidal inclusions embedded in a cube of matrix. 
The geometry of the problem is illustrated in Fig. [5] and details of the ellipsoids are 
given in Table [2| including the semi- axes' dimensions, centroid coordinates and Euler 
angles 0, u and ip^ which are successive rotations of the semi-axes ai, a2 and a^ about 
dobal coordinate axes z , X and z, respectively. The cube has side lengths of 600. 
The displacement boundary conditions were prescribed on faces x = 300, y = 300 
and z = 300 SlS Ux = 0^ Uy = and iZ^ = 0, respectively. The remaining faces at 
X = —300, y = —300 and z = —300 were subject to uniform normal unit tractions. 
The Young's modulus for the homogeneous matrix was chosen as = 1 and for the 
heterogeneities as = 2. Poisson's ratio was chosen as = 0.1 for both materials. 
All units are consistent. 

It is worthwhile noting that the small contrast in stiffness between the two ma- 
terials was chosen deliberately to maximise the extent of the perturbation ffelds 
emanating from the heterogeneities. Large contrasts in stiffness between the ma- 
trix and heterogeneities lead to perturbation fields that decay rapidly with distance 
from the heterogeneities. The close proximity of the three ellipsoidal heterogeneities 
to each other was also chosen deliberately in order to demonstrate the ability of 
the formulation to capture the interaction of multiple heterogeneities. Furthermore, 
the close proximity of one of the ellipsoids to a traction boundary was chosen to 
demonstrate the ability of the formulation to capture boundary effects. 



Incl. 


Centroid coordinates 


Semiaxes dim. 


Euler angles of max 


X 


y 


z 


ai 


a2 


0-3 








1 


-48.07 


78.27 


14.81 


50 


75 


100 


74.21 


48.44 


-48.07 


2 


16.45 


178.64 


-154.51 


50 


100 


75 


37.27 


22.27 


-25.51 


3 


127.93 


-65.94 


-27.32 


100 


75 


50 


46.74 


11.17 


-26.30 



Table 2: Topology and geometry of ellipsoidal inclusions of triple inclusion problem 



The problem was analysed using two three-dimensional finite element meshes 
with different densities, comprising C-HTS elements. The coarse mesh comprised 
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Figure 5: Geometry and topology of triple inclusion problem 



24 elements and 540 DOFs (Fig. [Gp) and the second, refined, mesh comprised 192 
elements and 3888 DOFs (Fig. lob). Results from the two enhanced finite element 




(a) (b) 

Figure 6: Triple inclusion task discretization by C-HTS elements: a) Coarse mesh with 24 enhanced 
elements (540 DOFs) b) Finer mesh with 1,536 enhanced elements (29,376 DOFs) 

analyses are plotted in the ^z-plane (at x = 0). Fig. [7^ and Fig. [8^ show the two 
meshes in this plane. The Gyy stress component for both analyses are shown in 
Fig. [7)3 and Fig. [8)3. 

In addition, a reference finite element analysis of the same problem was under- 
taken for comparison sake. The reference analysis utilised HTS elements but without 
the proposed enhancement. Unlike the other two analyses, the reference analysis 
utilised a mesh that explicitly resolved the three ellipsoidal heterogeneities and com- 
prised 309,406 tetrahedrons with 5,596,776 DOFs (Fig. [9^). The corresponding 
mesh and stress results of the reference analysis are shown in Fig. [9} Comparison of 
the stress results from the enhanced formulation and the reference analysis leads to 
the relative error plots shown in Fig. ^ and Fig. [8]3. It can be seen that even the 
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very coarse mesh with the enhanced formulation results in good agreement. 
Further comparison of the stress results is shown in Fig. 



10, where it can been 



seen that along the traction boundary, the finer mesh of enhanced elements is able 
to capture the imposed constant stress field more accurately than the very much 
finer reference mesh. The perturbation fields are based on the assumption of 




it 



0.010 

I 



Figure 7: Coarse mesh solution: a) Enhanced finite element mesh in yz plane at x = 0, b) Gyy in 
yz plane, c) error calculated as (cr^^^ — ^yyf' I 




a heterogeneity in an infinite medium but the enhanced formulation still exhibits 
convergence in the regions strongly infiuenced by the traction boundary. 

5.2. L-shaped specimen 

The proposed modelling strategy is also demonstrated on an example with a large 
number of inclusion. A 3D L-shaped specimen with fully fixed boundary conditions 
on the right surface of the right-hand arm and normal traction applied on the top 
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(a) (b) 




Figure 9: Reference analysis: a) discretization of entire body containing 309,406 tetrahedra, b) mesh 
refinement on surface of heterogeneities, c) Reference finite element mesh in yz plane at x = 0, d) 
crj:^^ in yz plane. 




(a) (b) 

Figure 10: Detailed 3I)-plots of cFyy stress concentrations due to the boundary effects in yz plane 
at X = 0: Comparison of solution from reference analysis with a) coarse mesh C-HTS solution with 
24 enhanced elements b) refined mesh C-HTS solution with 1,536 enhanced elements 
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surface of the left-hand arm is analysed, see Fig. [TT] The length of the plate is 300 in 
both X and y direction, the depth is 150 in z direction. The Young moduli were chosen 
as = 1 and E = 2 for matrix and inclusion respectively. Poisson's ratio was u = 0.1 
for both material phases. The microstructure comprised 2,523 spherical inclusions 
varying in size between 4 and 8 with a uniform spatial distribution (Fig. [TT})). All 
units are consistent. 

The solution for three different mesh densities (Fig. 11:, d & e) are compared 



in Fig. 12, These results are plotted on the x-y mid-plane. Fig. [12|(a-d) shows 
a plot of the axy stress component for the three different meshes. These results 
show that the complex stress distribution resulting from the heterogeneities can be 
captured and that the solution is converging with mesh refinement. Further local 
mesh refinement near corners and stress concentrations is possible, although this 
was not undertaken in this case. As an estimate of the computational overhead of 
the enhanced formulation for this particular problem on 16 processors, we note that 
the total solution for 30, 007 C-HTS elements was 597s, whereas the problem on the 
identical mesh of HTS elements for the equivalent homogeneous problem consumed 
1.6s of computer time. This represents a large increase of computational time in 
comparison to the homogenous problem. However, it should be noted that not all 
aspects of the solution procedure have been parallelised. It should also be noted 
that it was not possible to obtain the reference solution for this problem by means 
of conventional FEA due to the excessive number of degrees of freedom associated 
with a mesh required to resolve all of the heterogeneities. The mesh generation 
itself was not possible using our currently available software and hardware facilities. 
Furthermore, comparison with a PUM based solution, was also not undertaken due 
to the complexity of the problem. 



6. Conclusions 

A new micromechanics-enhanced finite element formulation has been presented 
for modelling the influence of a large number of heterogeneities in composite materials 
in a computationally efficient manner. The strategy exploits closed form solutions 
derived by Eshelby for ellipsoidal inclusions in order to determine the mechanical 
perturbation fields as a result of the underlying heterogeneities. Approximation 
functions for these perturbation fields are then incorporated into a finite element 
formulation to augment those of the macroscopic fields. A significant feature of 
this approach is that the finite element mesh does not explicitly resolve the hetero- 
geneities, although the resulting solution still explicitly accounts for their presence. 
In contrast with traditional homogenization approaches, this method does not rely 
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(a) (b) 




(c) (d) (e) 

Figure 11: L-shaped specimen a) geometry, b) microstructure comprising 2,523 spherical inclusions, 
c) coarse mesh comprising 1074 tetrahedral elements, d) medium mesh comprising 8772 elements 
and e) fine mesh comprising 30007 elements. 

on separation of scale and does not suffer from loss of information due to averaging 
or localization. 

The proposed technique has been implemented into a hybrid- Trefftz stress (HTS) 
element formulation and it has been shown that the resulting enhanced elements (C- 
HTS) require significantly fewer degrees of freedom to capture the detailed mechan- 
ical response compared to standard finite elements. The paper also outlines how the 
proposed micromechanics approach could be used within a Partition-of-Unity (PoU) 
formulation, although we conclude that this does not fully exploit the advantages of 
PoU methods and that the proposed hybrid- TreflFtz formulation is most appropriate. 

A self-compatibility algorithm is used to determine the mutual interactions be- 
tween inclusions, assuming that the eigenstrain fields are uniform within the domain 
of each inclusion. It was found that even for topologies exhibiting extremely small 
distances between the inclusions, this assumption is sufficient. Furthermore, it has 
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been shown that boundary effects, that are not accounted for by the classical mi- 
cromechanical solution due to the assumption of an infinite medium, can be captured 
through local mesh refinement. 
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We have implemented this formulation into our FE code that is optimized for 
parallel computing. Additional parallelization of the micromechanical aspects of 
the formulation needs to be investigated for increased efficiency. Further research 
is required in order to incorporate other improvements such as nonuniform eigen- 
strains [14j, debonding effects [24j and inclusions of arbitrary shape by evaluating 
the perturbation functions numerically [121 112 • 
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Appendix A. Detailed solution of perturbation displacements 

The displacement perturbation field in an infinite homogeneous material due to a 
uniform eigenstrain ejj applied to an ellipsoidal region Q is provided by the following 
integral equation [HI Eq. (11.30)] 



8^(1 - u) 



jkjki 



2v^kk,i - 4(1 - 



ik.kj 



(A.l) 



where u denotes the Poisson's ratio and the eUiptic potentials and $jj are defined 
as [H Eq. (11.32)] 



= ^[7 / |x - x'l dx' = ej.il), and = £[• / - — 3— dx' = 



(A.2) 



The integrals cj) and in Eq. (A.2) are the harmonic and bi-harmonic potentials 
respectively. Note that in Eq. (A.l) and thereafter, standard index notation is 
employed, together with the generalised summation convection due to Mura [14j. 
Thus, a repeated index is summed according to the Einstein summation rule (e.g. 
ciibij = Yl^=i^ibij)^ whereas a non-repeated upper-case index equals to the lower- 
case equivalent (e.g. aibiCjj = Yl^=i ^i^i^ij) • The symbol ajk^i denotes the partial 
derivative of ajk with respect to the coordinate x^. 
By combining Eq. (A.l) and Eq. (A.2), we obtain 



1 

8^(1 - u) 



(A.3) 
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Similarly to Eshelby's approach [Hj, the displacement perturbations are expressed 
in compact form: 



1 

87r(r^ jy) 



< = Ltjk^'jk. Lljk = ^ - 2z/5,,(/),, - 4(1 - iy)5,jct),k] (A.4) 



where the third-order operator Lijk maps a transformation eigenstrain 6:J^ to the dis- 
placement perturbation field . It is therefore analogous to the well-known Eshelby 
tensor [17j, which relates a transformation eigenstrain to the strain perturbation 
field. 

The operator Lij^ can be conveniently expressed in terms of the Ferrers-Dyson 
elliptic integrals, e.g. [14^ Eq. (11.36)] 

ds 



/(A) 
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4(A) 


27raia2as 



A{sy 

ds 



a] + s)A{sy 

(A.5) 



(af + s)(a2 + s)A(s) 
where stands for the i-th semi- axis of ellipsoid Q and A{s) is obtained from 

3 



A'{s) = ll{ai + sf (A.6) 



i=l 

The variable A is the largest positive root of equation [HI Eq. (11.37)] 

XiXi 



(ai + A)2 



= 1 (A.7) 



Notice that A is generally position dependent and non-zero for the points Xi placed 
outside the inclusion domain hence called the exterior points. Contrary, A = 
for interior points. 

All integrals in (A.5) admit a closed-form expression in terms of the Legendre- 
Jacobi integrals of the first and second kind, defined as a fuction of an auxiliary angle 
0, [H Eq. (12.17)]. It is worth noting that its definition via [14, Eq (11.18)] 



^ = sin-^/l-| (A. 
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is valid for interior points only and not everywhere as stated in [iMj. Thus, it needs 
to be replaced with a general formula: 



^ = sin- 



af + X 



(A.9) 



available, e.g. in [251E6]. Moreover, the following identity [IK Eq. (11.40.4)] 

[XnXnIi...jN{>^)]^p = 2XpIi,„jP + Ii...j^p{\) (A. 10) 

will be repeatedly proved useful in the sequel. 



It follows from Eq. (|A.4|) that to express operator L^^-^, we need to evaluate the 



first-order derivatives of the potential (j) and the third-order derivatives of To this 
end, we start with expressions [14^ Eq. (11.38)] 



0(A) = \ [/(A) - XnXjN{\)] 



and 



^,i(A) = \xi^2(t){X) - a] - x^x^//7v(A)]| = \xiQ{\) 
Employing Eq. ( |A.10D , the first derivative of (j) becomes 



(A.ll) 
(A.12) 



(t>,i{X) = \ |/i(A) - [x„x„/^(A)] ,} = \ [/,(A) - 2xih{\) - /,(A)] = -x,/,(A) 

(A.13) 



The third derivative of potential i\) is expressed from Eq. (A.12) as 



V^,^Jit(A) = \ \^^3QA^) + ^^kQA^) + (A.14) 

With the help of Eqs. ( A. IS] ) and ( |A.10 ), the term can be evaluated from 

g,,(A) = 2(/),,(A) - a\ [//(A) - XnXnhN{\)\^, = 2x, [ajlu{X) - /j(A)] (A.15) 
This provides the second derivatives of Q in the form 

QjkiX) = 2[Sjk [ajliAX) - Ij{X)] + Xj [ajluMX) ' Ij,kW] } (A.16) 

After utilising the derivatives of Q{X) and re-ordering the indices, Eq. ( |A.14 ) becomes 

^ifci(A) = XiSjk [aj/j/(A) - //(A)] + XkSji [a^IjKiX) - Ik{X)] 

+ xj5ki [a^jIjKiX) - Ik{X)] + XjXk [a]ljK,i{X) - IkA^)] (A.17) 
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with Iij.k and //j provided by ^ Eqs. (11.40.1, 11.40)] 

— 27raia2a3 



(a? + A)...(4 + A)(a2 + A)A(A) 



A 



2xp (gf + A)^ 



(A.18) 



Now we are in a position to evaluate Lf^-^. in terms of the Ferrers-Dyson integrals. 



Introducing Eqs. (A. 13) and (A. 17) into Eq. (A. 4) gives 



[87r(l - v)] Lf^.fc = XiSjk [a]lji{\) - //(A)] + {xkSji + XjSki) [a]ljK{\) - Ik{\)\ 

+ XjXk [a^IjKA^) - IkA^)] + 2v5jkXiIi{\) + 4(1 - z^)%a;fc/x(A)} 

(A.19) 



Since A = for the points inside the inclusion, recall Eq. (A. 7), all derivatives of /, 



vanish as well. Therefore, Eq. (|A.19|) yields 
[87r(l - v)] L^''"* 



XiSjk [a^jlji{\) - //(A)] + {xk5ji + Xj5ki) [aj/ji^(A) - /i^(A)] 
v5jkXiIi{\) + 4(1 - v)5ijXklK{\) (A.20) 



and Eq. (A. 4) receives its final form 



T£ _ T e.int I -*- 

ijk ijk ^87r(l-z^; 
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(A.21) 



For implementation purposes, it is worth noting that Eq. (A.4) admits the Voigt 
representation: 
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